Skyrmion Lattice in a Chiral Magnet 

S. Muhlbauer, 12 B. Binz, 3 F. Jonietz, 1 C. Pfleiderer, 1 * A. Rosch, 3 
A. Neubauer, 1 R. Georgii, 12 P. Boni, 1 

1 Physik Department E21, Technische Universitat Miinchen, D-85748 Garching, Germany 
2 Forschungsneutronenquelle Heinz Maier-Leibnitz (FRM II) 
Technische Universitat Miinchen, D-85748 Garching, Germany 
3 ITP, University of Cologne, Ziilpicher Str. 77, D-50937 Cologne, Germany 

*To whom correspondence should be addressed; E-mail: christian.pfleiderer@frm2.tum.de. 

February 11,2009 

Skyrmions represent topologically stable field configurations with particle-like 
properties. We used neutron scattering to observe the spontaneous formation 
of a two-dimensional lattice of skyrmion lines, a type of magnetic vortices, in 
the chiral itinerant-electron magnet MnSi. The skyrmion lattice stabilizes at 
the border between paramagnetism and long-range helimagnetic order per- 
pendicular to a small applied magnetic field regardless of the direction of the 
magnetic field relative to the atomic lattice. Our study experimentally es- 
tablishes magnetic materials lacking inversion symmetry as an arena for new 
forms of crystalline order composed of topologically stable spin states. 
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The formation of most crystals is related to three main aspects. First is the interplay of 
local repulsions and long range attraction of atoms leading to an instability of the liquid with 
correlations that are completely isotropic and without preferred direction. Second, three particle 
collisions lower the energy and dominate the formation of the crystal out of the isotropic density 
fluctuations. Third, atoms are quantized with an integer number of them in a unit cell. 

For the spin structures in magnetic materials it was believed that some or all of the three 
above mentioned mechanisms could not take place. Consider for instance the formation of 
a magnetically ordered state out of a paramagnetic metal. The underlying atomic lattice or 
the Fermi surface strongly breaks the rotational and translational invariance. Furthermore, for 
magnetic fluctuations in the paramagnet, three particle collisions are forbidden by time rever- 
sal symmetry. Lastly, the question arises, as to which quantized entities may play the role of 
the atoms in a solid. Possible candidates are topologically stable objects such as skyrmions, 
hedgehogs or merons, which have received great theoretical interest ranging from magnetic 
monopoles in particle physics (1, 2) to the emergence of gauge theories in condensed matter (3). 
In analogy to the Abrikosov vortex lattice in superconductors, these objects may be considered 
as the building blocks of novel forms of electronic order akin to crystal structures. However, 
the experimental evidence for these building blocks is scarce. 

We report the observation of the formation of a magnetic structure with hexagonal symme- 
try perpendicular to a small applied magnetic field in the cubic B20 compound MnSi. We show 
that this structure can be described approximately by a simple superposition of three helical 
states in the presences of a uniform field, where the superimposed state with the lowest energy 
can approximately be viewed as a lattice of antiskyrmion-lines, that is, magnetic vortices for 
which the magnetization in the center is antiparallel to the applied field. All three mechanisms 
mentioned above play an important role in explaining the magnetic structure. The lack of space 
inversion symmetry in the atomic crystal of MnSi results in weak spin-orbit coupling that gen- 
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erates slow rotations of all magnetic structures such that they decouple from the underlying 
atomic lattice very efficiently. Further, in the presence of an external magnetic field, which 
breaks time -reversal symmetry, three -particle interactions of the magnetic excitations do occur. 
Lastly, for the magnetic state that emerges under these conditions skyrmion lines, that is, certain 
topologically protected knots in the magnetic structure, take over the role of the atoms in usual 
crystals. 

At ambient pressure and zero applied magnetic field MnSi develops helical magnetic order 
below a critical temperature, T c = 29.5 K, that is the result of three hierarchical energy scales. 
The strongest scale is ferromagnetic exchange favoring a uniform spin polarization (spin align- 
ment). The lack of inversion symmetry of the cubic B20 crystal structure results in chiral spin- 
orbit interactions, which may be described by the rotationally invariant Dzyaloshinsky Moriya 
(DM) interaction. The ferromagnetic exchange together with the chiral spin-orbit coupling lead 
to a rotation of the spins with a periodicity \ h p=* 190 A that is large compared with the lattice 
constant, a ps 4.56 A. This large separation of length scales implies an efficient decoupling of 
the magnetic and atomic structures. Therefore, the alignment of the helical spin spiral along the 
cubic space diagonal (111) is weak and is only fourth power in the small spin-orbit coupling. 
These crystalline field interactions, which break the rotational symmetry, are by far the weakest 
scale in the system. 

Our study was partly inspired by recent work on the pressure dependence of the properties 
of MnSi (4, 5, 6). As shown in Fig. 1A well below T c an applied magnetic field, B, unpins the 
helical order and aligns its wave vector Q parallel to the field, Q || B, for a field exceeding 
B c i p^ 0.1 T (7, 8, 9). This state is referred to as the conical phase. For a magnetic field 
exceeding B c2 ~ 0.6 T, the effects of the DM interaction are suppressed, giving way to a spin- 
aligned (ferromagnetic) state. For temperatures just below T c an additional phase, referred to 
as the A phase, is stabilized in a finite field interval for B || (100) (Fig. 1) (7, 8). It had been 
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believed that the A phase was explained by a single-Q helix, where Q is perpendicular to the 
applied field (8). The specific heat exhibits a tiny peak at the border of the A phase, wheras the 
AC susceptibility discontinuously assumes a lower value (10). When taken together with the 
discontinuous change of the scattering intensity seen in neutron scattering, this suggests that the 
A phase represents a distinct magnetic phase, where the transition from the conical order to the 
A phase is discontinuous (first order). 

Because Q tends to align parallel to an applied magnetic field, neutron scattering as a func- 
tion of B has been reported for set-ups in which the magnetic field was perpendicular to the 
incident neutron beam (11). In contrast, we chose the incident neutron beam to be parallel to 
the applied magnetic field (Fig. IB). Two samples were studied. Sample 1 refers to a disk of 
19 mm diameter, d = 3 mm thick, where the vector normal to the disc was slightly misaligned 
by 11° with respect to a (110) axis. Sample 2 is a small parallelepiped, with dimensions 1.5 mm 
by 1.5 mm by 14 mm, where a (110) axis corresponded to the long axis. To search for higher 
order peaks and double scattering we increased the neutron flux in our measurements of sample, 
accepting a larger beam divergence. In contrast, the beam divergence was reduced for sample 
2 to improve the resolution in rocking scans. All data at finite magnetic field were measured 
after zero-field cooling to the desired temperature, followed by a field ramp to the desired field 
value. However, the results for the A phase were identical when recorded after field-cooling. 
For further details of the experimental set-up we refer to (11). 

Fig. 2 shows typical data recorded, where the spot sizes represent the resolution limit. All 
data shown represent sums over rocking scans of typically ±8° (11). Fig. 2A to C shows data for 
sample 1, whereas Fig. 2D to F shows data for sample 2. Fig. 2A shows the scattering intensity 
of sample 1 in a zero-field-cooled state at a temperature of 27 K for a (110) scattering plane, 
and the (110) axis is indicated in the figure. The pattern is consistent with previously published 
data and helical magnetic order along (111). Fig. 2B shows the intensity pattern of sample 1 in 
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the A phase. Six spots emerge on a regular hexagon. 

We tested the variation of intensity pattern on the orientation of the field relative to the 
crystal axes in both samples. The field was always parallel to the incident beam, whereas 
the sample was rotated for a large number of different orientations. Typical data are shown 
in Fig. 2C for sample 1, where the sample was rotated with respect to the vertical axis into a 
random position. For both samples and all crystal orientations the scattering pattern always 
exhibited the six-fold symmetry. In case the scattering plane contained a (110) direction, two 
of the peaks of the six-fold pattern coincided with this direction. As for Fig. 2C the scattering 
plane did not contain a (110) direction. For sample 2, the intensities along the vertical direction, 
which coincided with the (110) direction, were systematically weaker. This may be explained 
by the demagnetizing fields caused by the large aspect ratio, which implies that part of the 
scattering intensity was not captured in the rocking scans (see also (11)). The main result of our 
study is that for all orientations of the magnetic field with respect to the atomic lattice six Bragg 
reflections are observed on a regular hexagon that is strictly perpendicular to the magnetic field. 

We performed rocking scans to test whether the A phase has long-range order. Typical data 
are presented in (11). In the helical state the half-width of the rocking scans corresponded to 
a magnetic mosaicity rj m 3.5° consistent with previous work and long range order (12, 13). 
Remarkably, in the A phase the half-width of the rocking scans corresponded to a reduced 
magnetic mosaicity r\ m ps 1.75°, implying an even longer correlation length of at least £ ps 
5500 A, when allowing for demagnetizing fields (11). 

To test for consistency with previous work, we have also measured the emergence of the 
A phase as a function of temperature for magnetic field perpendicular to the neutron beam, 
where the vertical axis was the same (110) axis as before, and the low-symmetry horizontal 
axis containing spots 6 and 8 in Fig. 2F was perpendicular to the magnetic field and incident 
neutron beam. Data were recorded after: (i) zero-field-cooling the sample to a temperature well 
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below T c , (ii) increasing the magnetic field to 0.19 T and (iii) measuring the neutron scattering 
pattern for selected increasing temperatures (Fig. 2F shows data for T = 27.7 K). Well below 
T c we first observe the two spots parallel to the field direction labelled 9 and 10 characteristic 
of the conical state. When entering the A phase the intensity of the spots of the conical phase 
becomes very weak, but does not vanish, whereas strong scattering intensity appears in the 
perpendicular direction (spots 6 and 8). This is consistent with previous work and may signal a 
phase coexistence as expected of a weak first-order transition with possible extra effects of the 
demagnetizing fields added. 

The key results of our neutron scattering data may be summarized as follows: (i) the helical 
wave-vector aligns perpendicular to the applied magnetic field, (ii) the fundamental symmetry 
of the intensity pattern is sixfold suggesting a multi-Q structure, and (iii) the A phase stabilizes 
in a magnetic field strength of order B c2 /2. Moreover, the pattern aligns very weakly with 
respect to the (110) orientation. We can readily account for these features in the framework of 
standard Landau-Ginzburg theory by taking thermal fluctuations into account on top of a mean- 
field approximation. Near T c the Ginzburg-Landau energy functional can be written as (14, 15) 

F[M] = J d 3 r (r M 2 + J(VM) 2 + 2DM-(Vx M) + UM 4 - B • M) , (1) 

The first and second terms represent the usual quadratic contribution with the conventional 
gradient term, the third term the Dzyaloshinsky-Moriya interaction and the last term the cou- 
pling to an external magnetic field B. The quartic term accounts in lowest order for the ef- 
fects of mode-mode interactions and stabilizes the magnetic order. We neglect higher order 
spin-orbit coupling terms describing anisotropy effects (14, 15). The free energy is given by 
exp(— G) = /DM exp(— F[M]) (throughout the paper, we use a dimensionless free energy). 
Within mean field approximation G(B) is roughly equal to minF[M], and one minimizes F 
with respect to the spin structure M(r). 
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To explain the A phase, we evoke strong analogies with the crystal formation of ordinary 
solids out of the liquid state. The latter is in most cases driven by the cubic interactions of 
density waves (16), which in momentum space can be written as 

qi>q2,q3 

The ordered state can gain energy from this term only when three ordering vectors of the crystal 
structure add up to zero. Accordingly, in many cases (exceptions can arise only for strong first 
order transitions) the ordered phase which forms first out of a liquid state is of body-centered 
cubic symmetry (16), which is the crystal structure with the largest number of such triples of 
reciprocal lattice vectors. 

In the presence of a finite uniform component of the magnetization, Mj, a similar mech- 
anism can also occur in MnSi. From the quartic term in Eq. 1 we obtain terms which are 
effectively cubic in the modulated moment amplitudes 

( M /' m qi)( m q 2 •m q3 )5(qi + q 2 + q 3 ), (2) 

qi,q2,q3 

where m q is the Fourier transform of M(r). As in the case of an ordinary crystal, one can gain 
energy from this term for a structure with three Q-vectors adding up to zero. These vectors have 
a fixed modulus - determined by the interplay of the two gradient terms in Eq. 1 . Therefore these 
three vectors have relative angles of 120° (Fig. 3A) and define a plane characterized by a normal 
vector h. By symmetry, the energy change is proportional to Mj • h and therefore the three 
Q vectors must be perpendicular with respect to the external magnetic field. Our qualitative 
arguments already explain the two main experimental observations in the A phase: The Bragg 
spots are located in the plane perpendicular to B and display a sixfold symmetry (because both 
Q and — Q give a Bragg reflection) independent of the orientation of the underlying lattice. 
We therefore suggest that the A phase is a chiral spin crystal, the A crystal, approximately 
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characterized by the magnetization 
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(3) 



i=i 



where Mq.(r) = A(n n cos(Q;r) + n i2 sin(Qj)r) is the magnetization of a single chiral helix 
with amplitude A, wave vector Q; and two unit vectors, and n i2 , orthogonal to each other 
and to Qj. All three helices have the same chirality, i.e., all Qj ■ (n^i x n i2 ) have the same 
sign. More precisely, one has to add further higher-order Fourier components to Eq. 3 when 
minimizing F[M]. However, these terms remain small close to T c . The relative shifts Ar» of 
the helices, which we calculate theoretically, determine whether the A crystal can be described 
as a lattice of skyrmions (see below). 

One also has to take into account, that an external magnetic field favors helices with Q 
vectors parallel to B, because this way the spins may easily tilt parallel to the field to form a 
conical structure. Within mean-field theory, this conical phase always has the lowest energy 
(11). However, in the parameter range, where the A phase occurs experimentally, i.e., close 
to T c at intermediate magnetic fields (Fig. 1A), the energy difference between the two phases 
becomes very small as shown in Fig. 3B inset. The origin of the energy minimum of the A 
crystal for moderate magnetic field can be traced back to the size of the modulations of the 
magnetization amplitude, |M(r)| which is minimal close to B ps 0AB c2 (11)). In our mean- 
field Landau-Ginzburg theory, the A crystal thus appears as a meta-stable phase, which becomes 
extremely close in energy to the conical phase for intermediate fields B m 0AB c2 . 

It turns out that when we consider thermal fluctuations around the mean field solution these 
stabilize the A crystal. To show this, we consider the leading correction to mean field theory 
arising from Gaussian fluctuations 



where M is the mean field spin configuration for either the A phase or the conical phase. To 




(4) 
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make Eq. 4 well defined, one has to specify a cutoff scheme for short length scales. We use 
a cutoff in momentum space, k < 2-ir/a, where a is the lattice spacing of the MnSi crystal. 
Because the long pitch of the helix, it turns out that most contributions arise from fluctuations 
on short length scales with the exception of temperatures extremely close to T c (see (11) for a 
detailed discussion), but both short-range and long-range fluctuations favour the A crystal for 
intermediate magnetic fields. As shown in Fig. 3B inset, the fluctuations indeed stabilize the A 
crystal. A typical phase diagram resulting from (4) is shown in Fig. 3B. The theoretical phase 
diagram catches the main characteristics of the experimental phase diagram. The A crystal is 
stable at intermediate fields not too far from T c . When interpreting the theoretical result one 
has to take into account that Eq. 4 is only valid for small fluctuations and therefore can not 
be applied too close to T c . Indeed it is expected (17, 18) that fluctuations ultimately drive the 
transition first order and that such strong fluctuations will substantially shift the transition line 
to the paramagnet. We estimate the strength of fluctuations by calculating the leading correction 
to the order parameter for both the conical phase and the A crystal (11). In the shaded area of 
Fig. 3B these corrections are small (less than 20%), which justifies the use of Eq. 4. 

As can be seen from Fig. 3C, the magnetic structure of the A crystal obtained by minimizing 
F[M] is characterized by a pattern of magnetic vortices. To elucidate their nature we compute 
the skyrmion density given by (19): 

_ _ dn ^ dn ^ 

47T dx dy 

where x and y are the coordinates perpendicular to B and n = M(r)/|M(r)| is the orientation 
of the magnetization. is a measure of the winding of the magnetization profile. If integrates 
to 1 or —1, a topologically stable knot exists in the magnetization. As shown in Fig. 3D, 
the skyrmion density is finite and oscillates between positive and negative as compared with 
the normal helical or conical phases, where it is zero. Moreover, the skyrmion number $ = 
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/ 4>(r)d 2 r per two-dimensional unit cell is quantized and adds up to —1. Taken together this 
implies that the A crystal can be interpreted as a crystal made out of quantized objects, the 
skyrmion lines, with a magnetization at their core that is antiparallel to the applied magnetic 
field and M . 

Because of its symmetry and the cubic interactions, the A crystal has to be separated by two 
first order phase transitions both from the conical and the paramagnetic phases. This is consis- 
tent with the experimental observations. Two additional features in the scattering patterns that 
account for less than 1% of the total integrated scattering intensity are, first, weak continuous 
streaks of intensity emerging radially outward from the six main spots, and second, the coex- 
istence of conical order and the spin crystal. Both may be the result of weak heterogeneities 
resulting from these generic first-order boundaries of the A crystal; possibly in combination 
with demagnetizing fields. 

Many years ago, Bogdanov and collaborators used a mean-field model to establish the ex- 
istence of skyrmion lattices for anisotropic non-centrosymmetric magnetic materials under the 
application of a magnetic field (20, 21). The authors also pointed out that within their mean- 
field analysis for cubic materials such as MnSi the skyrmion lattice is always metastable. More- 
over, in the absence of a magnetic field, it has been shown that certain crystalline spin struc- 
tures can be stabilized by long-range interactions or an additional phenomenological parame- 
ter (22, 23, 24). In contrast to this work, we have reported here that it is sufficient to include 
the effects of Gaussian thermal fluctuations to stabilize skyrmion lattices in a magnetic field in 
cubic materials. 

It is instructive to search for analogies of the A crystal in other condensed matter systems. 
Because it is a multi-Q structure we note that previously known multi-Q structures, for exam- 
ple, in the rare earths involve large values of Q and exhibit very strong pinning to the atomic 
lattice (25, 26), whereas for MnSi we observe that the six-fold pattern of the A crystal exists in- 
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dependently of the underlying lattice and Q is quite tiny. Although flux lines in superconductors 
and the magnetic skyrmion lines observed are topologically completely different objects, there 
is nevertheless an intimate similarity of the Abrikosov lattice of superconducting flux lines and 
the hexagonal symmetry of the A crystal (20, 21). Moreover, the A crystal is characterized by 
broken translation symmetry in the plane perpendicular to B only. Therefore the A phase is 
similar to the chiral columnar phase of liquid crystals (16, 27). Further, the spin structure of the 
A crystal is topologically equivalent to theoretical predictions of the spin structure of the ferro- 
magnetic quantum Hall state near 1/2-filling (28), where, however, the underlying energetics is 
completely different. Lastly, individual magnetic vortices attract also great interest as a micro- 
magnetic phenomenon, which arises when conventional domain walls in soft ferromagnets are 
made to meet (29). 

The skyrmion lattice in the chiral magnet MnSi reported here represents an example where 
an electronic liquid forms a spin crystal made from topologically non-trivial entities. This 
provides a glimpse of the large variety of magnetic states that may be expected from the particle- 
like magnetic objects currently discussed in the literature (30). 
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Figure 1 : (A) Magnetic phase diagram of MnSi. For 5 = 0, helimagnetic order develops below 
T c = 29.5 K. Under magnetic field, the helical order unpins and aligns along the field above 
B c i, above B c2 , the helical modulation collapses. In the conical phase, the helix is aligned 
parallel to the magnetic field. The transition fields shown here have been inferred from the 
AC susceptibility, where the DC and AC field were parallel (100) (10). (B) Neutron scattering 
set-up used in our study; the applied magnetic field B was parallel to the incident neutron beam. 
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Figure 2: Typical neutron small angle scattering intensities; note that the color scale is loga- 
rithmic to make weak features visible. Data represent the sum over rocking scans with respect 
to the vertical axis through the sample. (A) to (C) show data for sample 1 and (D) to (F) for 
sample 2. Backgrounds measured above T c and B = were subtracted in all panels except 
for (A) (light blue square). Spots are labelled for reference; for the intensity of these spots as 
a function rocking angle see Ref. (11). (A) Helical order in sample 1 in the zero-field-cooled 
state at T = 27 K and B = 0; (B) Six-fold intensity pattern in the A phase in sample 1; same 
orientation as in (A); T = 26.45 K, B = 0.164 T (C) Six-fold intensity pattern in the A phase 
for random orientation of sample 1 (see text for details); T = 26.77 K, B = 0.164 T (D) He- 
lical order in sample 2 in the zero field cooled state at T = 16 K and B = 0; (E) A phase in 
sample 2, same orientation as in (D); T = 27.7 K, B = 0.162 T (F) A phase as measured in 
conventional set-up (compare Fig. 1A in (11)), where data in all other panels were measured in 
the configuration shown in Fig. IB); T = 27.7 K, B = 0.190T. A small residual intensity due 
to the conical phase is observed (spots 9 and 10), whereas spots 6 and 8 correspond to those in 
(E). 
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Figure 3: (A) Depiction of the hexagonal basis vectors of the crystalline spin order in the 
A phase. (B) Theoretical phase diagram as a function of magnetic field B / B with B Q = 
J (JQ 2 ) 3 /U and the parameter t = r Q J/D 2 — 1, which is roughly proportional to T — T c . 
We use the model parameter 7 = JD/U = 5 and a momentum- space cutoff k < AOD/J. 
Smaller values of 7 increase the A phase regime. For most values of field, the A phase is 
either metastable or stable, but at low fields below the dotted line, it becomes unstable. Above 
and to the right of the red dashed line, fluctuation correction to the size of the order parameter 
becomes larger than 20% and our theoretical analysis becomes uncontrolled. Therefore in the 
shaded gray region, we have therefore reliably established stability of the A phase within our 
model (see (11) for details). (Inset) Energy difference between A phase and conical phase as a 
function of field for the same parameters and t = —3.5, both in the mean-field approximation 
and with fluctuation corrections. Fluctuations stabilize the A phase at intermediate fields. (C) 
Real space depiction of the spin arrangement in the A phase in the x-y plane. Note that this 
spin arrangement is translation-invariant along the z-axis, which is parallel to the magnetic 
field. (D) Skyrmion density per unit cell area as calculated for the A phase as shown in panel 
(C). The integrated skyrmion density per unit cell is finite, $ = — 1. The arrows represent the 
magnetization component perpendicular to the line of sight. 
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1 Experimental details 



Our neutron scattering measurements were performed at the diffractometer MIRA at FRM II at 
the Technische Universitat Miinchen (SI). Data were recorded for an incident neutron wave- 
length A = 9.6 A with a 5% FWHM wavelength spread. A delayline 3 He area detector of 
200 x 200 mm 2 was used with a position resolution of order 2 x 2 mm 2 . The neutron beam 
was collimated over a distance of 1.5 m. For sample 1 the source aperture and the aperture at 
the sample were both 8x8 mm 2 . For sample 2 the source aperture and the aperture at the sam- 
ple were both 4x4 mm 2 . However, the size of sample 2 of 1.5 x 4 mm 2 acted as an effective 
smaller aperture at the sample reducing the beam divergence even further. The distance between 
the sample and the detector was between 0.8 and 1.3 m. 

Samples were cooled with a cryogen free pulse tube cooler. The magnetic field was gener- 
ated with bespoke water-cooled Cu solenoids in a Helmholtz configuration (52). The magnetic 
field profile was carefully characterized with a Hall probe and found to be uniform much better 
than 1% over the sample volume. All data at finite magnetic field were measured after zero-field 
cooling to the desired temperature, followed by a field ramp to the desired field value. However, 
in the A-phase data were identical also after field-cooling. 

Because the ordering wave-vector Q in MnSi over large portions of the magnetic phase 
diagram tends to align parallel to an applied magnetic field, neutron scattering as a function 
of B has been reported for set-ups where the magnetic field was perpendicular to the incident 
neutron beam as shown in Fig. SI A. In contrast, we chose the incident neutron beam to be 
parallel to the applied magnetic field as shown in Fig. SIB. 

Two samples were studied. Sample 1 refers to a disk of 19 mm diameter, d = 3 mm thick, 
where the vector normal to the disc was slightly misaligned by 11° with respect to a (110) axis. 
Sample 2 refers to a small parallepiped, with dimensions 1.5 mm x 1.5 mm x 14 mm, where a 
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(110) axis corresponded to the long axis. 

2 Demagnetizing fields 

In our study we find that the hexagonal magnetic scattering intensity in the A-phase aligns 
strictly perpendicular to the applied magnetic field. This implies that the magnetic scattering 
intensity depends sensitively on the distribution of magnetic field directions across the sample 
volume, notably the effects of demagnetizing fields. 

The perhaps most pronounced consequence of the demagnetizing fields is seen in the ori- 
entation of the hexagonal scattering intensity with respect to the orientation of sample 1 . As 
stated above this sample is a circular disc of 19 mm diameter, d = 3 mm thick, where the vector 
normal to the disc was slightly misaligned by 11° with respect to a (110) axis. When the ap- 
plied field was not perfectly perpendicular to the disc, we observed a deflection of the maxima 
consistent with the demagnetizing fields as shown in Fig. S2. 

3 Rocking scans 

In our study rocking scans could only be carried out for the vertical axis. This is particularly 
important for sample 2, which was a bar with a cross-section of 1.5 x 1.5 mm 2 and 14 mm 
long. The long side of sample 2 was parallel to a (110) direction. In the neutron scattering 
experiment the long side of the sample was vertical and perpendicular to the incident neutron 
beam and magnetic field (cf set-up in Fig. S IB). Thus in the A-phase two spots of the hexagonal 
scattering intensity coincided always with the vertical (110) axis and thus the axis of the rocking 
scans (cf Fig. 2E in the main text). 

We note that for a reasonably well collimated neutron beam with a divergence smaller than 
the rocking width, the observation of equal intensity of the vertical spots implies, that neither of 
these vertical spots satisfies the Bragg condition. The intensity corresponds only to the tails of 
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the Bragg spot. Therefore, when summing over rocking scans with respect to the vertical axis, 
these vertical spots remain weak. In other words the scattering intensity in the vertical spots was 
not fully captured in rocking scans, with possible small effects of demagnetizing fields added. 

We have also observed weak intensity at the positions of higher-order peaks, although this 
may arise from double scattering rather than true higher order reflections. An example of double 
scattering is the appearance of weak spots along (100) in Fig. 2A of the main text as expected 
of a multi-domain state of helical order at B = 0. As shown Fig. 2B and C of the main text to 
achieve a higher flux and thus larger intensity in sample 1, which had a larger sample volume, 
we accepted a larger beam divergence. This revealed additional weak spots for the A-phase. 
However, because the overall intensity was still very weak, it was not possible to distinguish 
unambiguously if these spots were due to double scattering alone or whether they represented 
higher order peaks. 

While the search for higher order peaks was inconclusive, it was possible to infer evidence 
of long range order from the rocking scans shown in Fig. S3. In the helical state the half- width 
of the rocking scans corresponded to a magnetic mosaicity r\ m ss 3.5° consistent with previous 
work and long range order (S3). Remarkably, in the A-phase the half-width of the rocking scans 
corresponded to a reduced magnetic mosaicity r\ m rs 1.75°, implying an even longer correlation 
length of at least £ ps 5500 A when allowing for demagnetizing fields, where the Lorentz factor 
was taken into account. A more refined analysis of the peak shape, e.g., considering issues 
discussed in Ref. (S4), is beyond the present study. 
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4 Theoretical model 



In the following we give a detailed description of our theoretical calculations where we show 
that thermal fluctuations can stabilize the skyrmion crystal. We begin with a few remarks of pre- 
vious work on skyrmion lattices in Section 4.1. Due to the long and rich history of the subject in 
fields ranging from high-energy physics to soft matter the discussion will unfortunately remain 
incomplete. We next revisit the standard phenomenological Landau-Ginzburg model for heli- 
magnets such as MnSi in an applied magnetic field. We analyze this model taking into account 
Gaussian fluctuations on top of the mean-field solution. Model and approximation scheme are 
introduced in Section 4.2. The well-known mean-field ground state of this model is the conical 
phase. In Section 4.3, we give simple arguments that for intermediate fields, there is a compet- 
ing ground state, which corresponds to a triangular spin crystal of skyrmions (the A-crystal). 
In Section 4.4, we consider the mean-field approximation and find that while the conical phase 
remains the mean-field ground state, the energy of the A-crystal comes very close at intermedi- 
ate fields. We then show in Section 4.5 that the fluctuation corrections to the free energy lower 
the energy of the A-crystal compared to the conical state and make it stable. Finally, in Section 
4.7, we comment on the experimentally observed very weak pinning of the magnetic A-crystal 
to the atomic crystal. 

4.1 General remarks 

Continuous fields with topologically non-trivial, stable, particle-like properties have attracted 
great interest in a large number of disciplines for many decades (see also (55)). Many of these 
studies have been inspired by the Skyrme's seminal work (S6). He was able to show that in 
a certain field theory describing interacting mesons one can find quantized and topologically 
stable field configuration. Surprisingly these excitations made of bosonic fields can be inter- 
preted as baryons, i.e. as fermions. In the area of nuclear and particle physics, a very large body 
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of work has been published related to Skyrme's work. Regarding the formation of skyrmion 
lattices and related structures it has, for instance, been suggested that dense nuclear matter in 
neutron stars may be described by a three-dimensional lattice of skyrmions (57), where typical 
face-centered cubic structures give the lowest energy (S8, S9). 

Skyrmion lattices and related structures have also been considered in many other areas of 
physics. A long time ago, Bogdanov and coworkers (S10) showed for non-centrosymmetric 
magnetic materials that chiral interactions represent an elegant route to skyrmion ground states 
(S10, Sll). In these studies it was found theoretically that for a class of anisotropic systems 
skyrmion lattices are the mean-field ground states for a range of parameters. Experimental ex- 
amples confirming these predictions remain to be discovered. The authors also pointed out the 
similarity of these states to vortex lattices in superconductors. Indeed, it was emphasized by 
Bogdanov and Hubert (S12), that a magnetic field can be instrumental in stabilizing skyrmion 
lattices. The authors also pointed out that, within their mean-field analysis for cubic materials 
like MnSi, the energy of a conical state is always lower than the one of skyrmion lattices as 
is discussed in Section 4.4 below. One way to stabilize crystalline magnetic structures, e.g., 
the skyrmion lattice, for cubic systems on the mean-field level is to consider long-ranged in- 
teractions (S13, S14) or extra phenomenological parameters (S15) added to the conventional 
Ginzburg-Landau approach. Instead, we will show below, that it is sufficient to include the 
effects of Gaussian thermal fluctuations to stabilize skyrmion lattices in a magnetic field. 

A large body of literature on skyrmion lattices is also available in the area of quantum 
Hall systems. For instance, the two-dimensional analog of the Skyrme model can be used 
to describe excitations in ferromagnetic quantum Hall systems with small Zeemann splitting 
where the spin-density takes over the role of the pion field (S16). In these systems the quantized 
winding number of the skyrmion can be identified with the charge of the excitation. Several 
theoretical proposals suggested also lattice states of these two-dimensional skyrmions (S17,S18, 



6 



S19). Nuclear magnetic resonance measurements (S20, S21) and recent inelastic light scattering 
experiments (S22, S23) gave results consistent with the existence of such a skyrmion lattice in 
GaAs heterostructures with fillings close to v — 1, but there is no direct microscopic evidence. 

The physics of magnets without inversion symmetry finally shares many similarities with 
cholesteric liquid crystals. In these cholesterics, a sequence of several phase transitions from 
helical to so-called 'blue phases' has been observed. Blue phases are characterized by complex 
patterns of order parameters woven from topological defects of the underlying chiral helices 
(S24) which appear as colorful objects when typical lattice distances are of the order of the wave 
length of visible light. Skyrmion textures have also been considered in these systems (S25). 

4.2 Definition of the model and saddle point approximation 

It has long been established (S26, S27), that the helimagnetism in MnSi is well described by 
the following leading-order Landau-Ginzburg functional, which depends on the continuously 
varying magnetization M: 



F[M] = / d 3 r (r M 2 + J(VM) 2 + 2D M • (V x M) + UM 4 - B • M) , (SI) 



where B is the external magnetic field and r ,J,D,U are parameters (U, J > 0). Here, we 
chose D > 0, which selects a left-handed spiral with wavevector Q = |Q| = D/J. Within 
the Landau-Ginzburg approach, one linearizes all T dependences around T c . We therefore keep 
only a linear T-dependence of r (see below). 

The (dimensionless) free energy Gas a function of magnetic field and temperature is ob- 
tained as 



To evaluate the functional integral in Eq. (S2), we use the saddle point approximation (method 
of the steepest descent, Laplace's method), which consists of expanding F around its local 
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(S2) 



minima and performing Gaussian integrals over fluctuations around these minima (S5). For a 
given local minimum M (r), the result is 



where we have omitted an additive constant. The first term of the right-hand side is the mean- 
field contribution G mf and the second term, GWt, is the correction from Gaussian fluctuations. 
If corrections from fluctuations become of order 1 in the Ginzburg regime close to T c , the saddle 
point approximation breaks down and Eq. (S3) is no longer valid. 

By choosing appropriate units, we can eliminate two of the four parameters in Eq. (SI). For 
this, we re-scale all lengths as f = Qr, magnetizations as M = [U/( JQ 2 )} 1 ^ 2 !^, and fields as 
B = [U/iJQ 2 ) 3 } 1 ^ and obtain 



wheret = r /(JQ 2 ) — 1 oc T— T C MF measures the distance to the mean-field critical temperature 
T C MF (i.e. within the saddle-point approximation, the system is spiral spin-ordered for t < 
and paramagnetic for t > 0) and 7 = J 2 Q/U provides a relative weight between the mean-field 
and fluctuation contributions to Eq. (S3). Therefore, the physics of the model (SI) depends on 
three parameters: 7, t and the magnetic field. From now on, we will omit the tildes in most 
formulas to simplify notation, but keep in mind that we have chosen particular units. 

4.3 Intuitive arguments for the skyrmion lattice solution 

It is possible to give simple geometric arguments, which lead to the specific spin crystal adver- 
tised in this paper. We may single out the ferromagnetic component of the magnetic structure 
M/ = / M(r)d 3 r/V, where V is the volume. The ferromagnetic component M/||B is induced 
by the external field. The true order parameter is therefore $ = M — Mj. The quadratic part 




(S3) 



F = 7 \Sr [(1 + t)M 2 + (VM) 2 + 2 M • (V x M) + M 4 - B • Ml , 




(S4) 
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is minimized by a helix with wave length 2-k/Q described by 

M^(r) = A(ni cos[Qr] + n 2 sin[Qr)], (S5) 

where the wave vector Q and two unit vectors, r^, n 2 are orthogonal to each other. The sign of 
D in Eq. SI determines the chirality with Q • (ni x n 2 ) > for D > 0. To quadratic order, also 
arbitrary linear combinations of such helices minimize the free energy and only the interaction 
M 4 selects the magnetic structure. Expanding this term, we obtain 

M 4 = M) + 4Mj $ M f + 2Mj$ 2 + 4{M f ■ $) 2 + 4$ 2 $ • M f + $ 4 (S6) 

The crystaline state can gain energy from the second last term oc $ 2 $ ■ M/, which is cubic in 
the order parameter. Fourier transformation leads to 

|M / -$$ 2 rf 3 r= Yl (M / .m qi )(m q2 -m q3 )5(q 1 + q 2 + q3), (S7) 
qi,q2,q3^o 

where m q is the Fourier transform of M(r). Therefore, the cubic term vanishes unless the 
magnetic structure contains Fourier modes of at least three wavevectors Qj with Qi+Q 2 +Q3 = 
0. If in addition, |Q^| = Q, we immediately arrive at the coplanar arrangement of Fig. 3 A of 
the main text. Thus we obtain six Bragg spots with 60 degrees angles between them. Direct 
calculation shows that in this situation, the vector / $ 2 $ is orthogonal to the plane containing 
Qi, Q 2 , Q 3 with a prefactor that depends on the relative phases of the three helices. Hence, 
the cubic term is minimized if Q 1; Q 2 , Q 3 are orthogonal to B (the relative phases then ensure 
that / $ 2 $ becomes anti -parallel to M/, which minimizes the energy gain by the cubic term). 
Hence, minimization of the cubic term already explains the experimentally observed Bragg 
spots with hexagonal symmetry in the plane orthogonal to the field. 

We have also checked quantitatively that the cubic term plays an essential role in the ener- 
getics of the crystal phase. Namely, if the cubic term is artificially put to zero in the evaluation 
of the free energy, the crystal state no longer comes energetically close to the conical phase at 
intermediate fields. 
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4.4 Calculation of the mean-field free energy in an external field 

The mean-field ground state (i.e. the global minimum of the functional F[M]) can be found 
rigorously. For B < B c2 = (— 2t) 1 / 2 , it is the well-known (S28) conical phase M(r) = mz + 
a(x cos(z) + y sin(z)) with m = B/2, a = \\f— 2t — B 2 and with the free energy 

^° f n = -^(t 2 + B 2 ), (S8) 

where V is the system volume. The conical phase evolves continuously into the helical state at 
zero field. For larger fields, B > B c2 , the ferromagnetic state becomes the ground state. The 
proof is provided most easily by writing F[M] in the following way: 



lF = -\/^^-+\/E m - q [^ /3 (q)-^ /3 ] m q+/( M2 + d3r + v [ M f-*j » 

^° (S9) 
where r ab (k) = (1 + t + k 2 )5 ab — 2ie abc k c . The conical phase manages to minimize each term 
of the right-hand side individually. 

However, as argued above, there is a competing state at intermediate fields, which we call 
the A-crystal. This state can be easily addressed in the limit of small t < 0, where it is obtained 
from a superposition of the uniform component M/ with three left-handed helical spin spirals. 
More precisely, 

3 

M(r) = M f + (m Q ,e^ r + c. c.) , (S10) 

where Qi, Q 2 , Q3 are wavevectors orthogonal to the field, with |Qj| = Q and mutual angles of 
120° between them as shown in Fig. 3A of the main text, "c. c." denotes the complex conjugate 
and 

1 

m Q 3 = 2^ ^Qj + ^Qj ) ' < S 1 

where e'q., €q. are orthogonal unit vectors with x €q. = Qj and ipj is a complex number 
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which encodes the amplitude and phase of the helical spin-density wave j. Thus, 

M(r) = M f + £ |Vil(cos(Q,- • r + <^)e' Qj - sin(Q, • r + ^)e'^), (S12) 

where ^ = |^-| exp (i(f)j). Higher order Fourier modes vanish in the limit of small —t (but 
are included in our numerical calculations for finite t). The minimization process leads to 

I = 1^2 1 = 1^3 ]■> i- e - me three helices have equal weight. In addition, minimization also 
fixes the relative phases of the three helices in such a way that at one point (say r = 0), the 
magnetization of each helix points opposite to the field direction. This point is the center of 
the anti-skyrmions, shown in Fig. 3C of the main text. The values of the relative phases have a 
strong influence on the resulting magnetic structure and determine whether one obtains a lattice 
of skyrmions as described in the main text. 

In Fig. S4, we plot 

^AG mf =^(G mf -GS?) (S13) 

as a function of field in the limit of small t < 0. The energy difference to the conical phase 
is therefore plotted in units of 7t 2 /4 which corresponds to the energy difference between 
ferromagnetic and conical states at zero field. For comparison, we have also included (i) 
the energy of a single helix oriented orthogonal to the field, which corresponds to the tra- 
ditional interpretation of the A-phase, and (ii) the energy of the ferromagnetic state. The 
free energy difference between the A-crystal and the conical state obtains its minimal value 
(4/7tV)AG = 0.02838 at an external field of B = 0.403575^. At this point, the A-crystal 
obeys M f = 0.7773 1 = B/2. 

As can be seen in Fig. S5, the field where the energy of the A-crystal and the conical phase 
are closest, is characterized by small spatial variations of the magnetization amplitude |M| . This 
observation reflects the fact that it is energetically unfavourable to have too strong modulations 
of the magnetization amplitude. In particular, it would cost a large energy to suppress M locally 
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to zero, which makes it meaningful to classify magnetic configurations into topological sectors. 
For example, we find that the A-crystal and the conical phase belong to different topological 
sectors, since both configurations cannot be continuously transformed into each other without 
suppressing M locally to zero. In contrast, other non-trivial spin textures which have been 
proposed earlier in the context of noncentro symmetric magnets (S29, S15, S14), all have points 
or lines with vanishing magnetization, making any notion of topological stability ill-defined. 

Equation (S12) gives the exact form of the A-crystal in the limit t — > 0, very close to the 
critical temperature. But for t < 0, the crystal obtains additional corrections, which further 
suppress amplitude modulations. Namely, the spin structure obtains Fourier weight from all 
wavevectors of the reciprocal lattice nQi + mQ 2 and each wavevector is composed of two 
transversal modes (left- and right-handed helices) and a longitudinal one (where rnq || Q). To 
calculate the mean-field energy, we include wavevectors up to a short-distance cutoff, |Q| < A, 
and minimize numerically with respect to amplitudes and phases of all modes. The structure of 
the A-crystal gets therefore distorted from the simple form of Eq. (S12), but without changing 
its symmetry or topology. 

In our mean-field Ginzburg-Landau theory, the A-crystal thus appears as a metastable phase, 
which for intermediate fields becomes remarkably close in energy to the conical phase. Exper- 
iment suggests that close to the critical temperature, the crystal state becomes stable for these 
fields. In the following Section, we show that fluctuation corrections to the mean-field theory 
readily explain the stability of the A-crystal within the model (SI). 

4.5 Fluctuation contribution to the free energy 

We now consider the effect of fluctuation contributions to Eq. (S3): 




(S14) 
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where g is the matrix of second derivatives, 

1 d 2 F 

9 ™ = 2dm\dml, (S15) 

= 7(5k,k'r a6 (k) + 25 a6 ]Tm^ (S16) 

\ k" k" / 

We will be interested in the energy difference between the conical phase and the A-crystal 

AGfl uc t = Gfl uct — Gfl°cf (S17) 

Surprisingly, the energy difference between conical phase and A-crystal obtains a contribu- 
tion even from short length- scales as we show in the following. First, we write g^, = g + g±, 
where g is the first term on the right-hand side of Eq. (S 16), which is independent of the mag- 
netic state, and g\ are the two remaining terms, which are quadratic in m k . For large values of 
k or k', the term g behaves asymptotically like g m 5 Q/3 5kk' k 2 and dominates over g 1 . The 
short- wavelength contribution can therefore be expanded as 

G fluc t = - (trlog^ + tr^" 1 ^) + tr(g 1 g 1 g 1 g 1 ) + ...). (S18) 



We finally obtain the short- wavelength contribution to the energy difference as 

(S19) 



I 

_ A ^short _ UJV 

\/ AGfluct _ 2n 2 



- + - I M- 
2 V 



where A is an ultraviolet wavevector cutoff, M j4 (r) is the magnetization in the A-crystal and 
we have used the fact that M 2 = — 1/2 in the conical phase. In the presence of higher-order 
gradient terms in g one can use the equation / g d 3 k/(2n) 3 = A/(2tc 2 ) as a definition of A. 
The remaining part 

AG l °2 = AG fluct - AGS (S20) 

has a finite limit for A — > oo and needs to be calculated numerically. For periodic structures 
such as the conical phase or the A-crystal, one may decompose k = Q + q, where Q is a 

13 



reciprocal lattice vector and q is in the first Brillouin zone. The fluctuation matrix is then of the 
form g&, = <J q>q , g$ Q ,(q), such that 



where the integral is over the Brillouin zone. In principle, we could chose a one-dimensional 
reciprocal lattice for the conical phase and a two-dimensional one for the A-crystal. However, in 
order to treat the two phases on equal footings, we chose a three-dimensional reciprocal lattice 
spanned by Q 1; Q 2 from Fig. 3A in the main text perpendicular to the field and a third vector 
of the same length along the field. Thus, the Brillouin zone obtains the shape of a hexagonal 
prism. We introduce the ultraviolet cutoff by truncating the matrix at |Q|, |Q'| < A. 

To obtain the correct Gaussian fluctuations it is important to use the numerically exact 
saddle-point within the chosen cutoff. Therefore, we can not restrict it to the simple form of 
Eq. (S 12) for the A-crystal but perform a numerical minimization with respect to all modes rriQ 
for | Q | < A to determine the saddle point. We also minimize with respect to the fundamental 
helix wavevector length Q, which turns out to be slightly reduced from unity. 

After the minimization, which determines AG m f, we calculate / M 2 and obtain AG^/. 
Then, we calculate the fluctuation matrices and evaluate the Brillouin zone integral numerically, 
to determine AGjj^. In Fig. S6, we show, how the different contributions to AG = AG mf + 
AGj h u ° c f + AG^™! depend on the cutoff. These results confirm that AGjj° u n c g t has indeed a finite 
limit for A — > oo. These figures also show that fluctuation contributions at both short and long 
length- scales favor the A-crystal over the conical phase. Due to the thermal fluctuations the 
A-crystal obtains energetic stability as shown in the inset of Fig. 3B in the main text, where 
AG m f and AG are shown in units of 7t 2 /4. 

In real systems, the cutoff is of the order of 2n/a, where a is the lattice constant. For the 
case of MnSi, we have A AOQ (where Q = 1 in our units). For such large cutoff values, the 




(S21) 
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A ^ A ^ 5A 

AG « AG mf + — 



(S22) 



short-distance contributions dominate over the long-distance ones, such that 

is a very good approximation except for extremely small values of \t\ (the asymptotic small-t 
behavior is AG mt - oc t 2 , AG S ^ oc t and AGfl u n | oc This approximation was used to 

establish the phase diagram shown in Fig. 3B in the main text. 

Interestingly, the A-crystal is even no longer a local minimum of F for too small values of 
B, i.e., below the red dotted line in Fig. 3B in the main text. There, it becomes locally unstable 
(i.e. unstable with respect to small variations) since the fluctuation matrix 5 2 F/(5~M) 2 obtains 
a negative eigenvalue. However, we find that the A-crystal at intermediate fields is both locally 
stable and energetically favorable to the conical phase due to fluctuation corrections to the free 
energy. 

4.6 Testing the validity of the saddle-point approximation 

The saddle-point approximation, as expressed in Eq. (S3), is obtained by an expansion around 
a local minimum of F. It is thus only valid, if fluctuations around the minimum are small or, 
technically, if fluctuation corrections to the size of the order parameter are small compared to its 
mean-field value. To test this hypothesis, we calculate the response with respect to a fictitious 
space-dependent field B(r) and obtain the local magnetization as M(r) = — 5G/SH(r). Using 
Eq. (S3), we obtain 

M ' = M »-k n SMWM^ X '''\ m ' <S23) 

JVLo 

where M l = M a (r) (i.e. i, j indices stand for both a space coordinate as well as a directional 
index) and x 1 — l/25 2 F/(5M) 2 . The second term on the right-hand side is the fluctuation- 
correction to the order parameter. If the saddle-point approximation is valid, we expect this 
correction to be small. As a test for the validity of our approximation, we therefore demand that 
the second term in Eq. (S23) is small compared to the first one. 
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More precisely, we take the Fourier transform of Eq. (S23) and obtain the fluctuation- 
correction to the order parameter as 

Qi,Q2,Q3 ' ^ ' 

(S24) 

where Xqq^q) is the inverse matrix of <?QQ/(q). For the saddle-point approach to be valid, we 
demand that 

|<Jm Q | < |m Q | (S25) 

for the dominant modes with | Q | = 1 . 

In the same way as above for the fluctuation-correction to the free energy, it turns out that 
Eq. S24 obtains a cutoff-dependent contribution from short length-scales as follows: 

SA _ 

5m Q = -~ZaZ E ™-Q'XWO) + 0(A°). (S26) 

Z7F ' Q',6 

For large values of the cutoff, such as A 40, the short-distance part (oc A) dominates over 
the long-distance part [0(A°)]. We therefore use Eq. (S26) to map out the regime of validity, 
which we define as the region, where |5mq| < 0.2|mq|, in the phase diagram [red dashed line 
in Fig. 3B in the main text]. 

4.7 Pinning of spin crystal to the atomic lattice 

In our experiment we observe that the six-fold diffraction pattern of the A-phase gets weakly 
oriented with respect to the atomic crystal within the plane perpendicular to the magnetic field. 
The pinning to the crystal is such that two of the six peaks point in a (110) crystal direction, as 
long as this is compatible with the field orientation. Here we show that this behavior is fully 
compatible with our phenomenological theory of the A-crystal. 

The model of Eq. (SI) is fully rotation symmetric. A rotation of the magnetic structure is 
performed by M(r) — > KM(R^ 1 r), where R is an 5*0(3) operator. The coupling to an external 
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field breaks the rotation symmetry to SO (2), the group of rotations around the field-direction. 
Terms which are higher-order in spin-orbit coupling (neglected in Eq. (SI)) will break this 
symmetry and couple the magnetic structure to the underlying atomic lattice. 

To investigate pinning systematically, we organize these terms in powers of spin-orbit cou- 
pling (S30). In our mean-field Ginzburg-Landau approach, this means t — > and M ~ B* ~ 
(— t) l l 2 . The energy scale of helix formation is of second order in spin-orbit coupling (S26,S27). 
The leading anisotropy terms are of fourth order and given by (S31) 

= Ai £(g 4 + t v + 9 2 4 )|m q | 2 + A 2 / M 4 + M 4 + M 4 . (S27) 

q J 

We treat these anisotropy terms perturbatively, i.e., we first determine the magnetic structure 
using Eq. (SI) and then calculate how AF^ 4 ^ changes if the structure is rotated around the field 
direction. This gives a pinning function F(ip), where if is the rotation angle relative to some 
standard orientation. 

The leading anisotropy term of Eq. (S27) provides the pinning of the single helix state in 
MnSi and related helimagnets. In MnSi, the parameters Ai, A 2 are such that g 4 + g 4 + g 4 is 
minimized, i.e., the preferred spiral direction at zero field is (111). If we assume for a moment, 
that the A-phase consists of a single spiral which is oriented orthogonal to the field, this spiral 
would lock into a direction q which minimizes g 4 + g 4 + g 4 under the constraint that q • B = 0. 
This scenario is not compatible with the observations from neutron scattering (Fig. 2 in the main 
text). 

Interestingly, the anisotropy term in Eq. (S27) does not lead to a pinning of the hexagonal 
spin crystal to linear order in Ai and A 2 (i.e. to fourth order in spin-orbit coupling). The reason 
is the six-fold rotation axis of this crystal. As a consequence, the pinning function is of the form 
F{v) — J2n( a n cos 6n(p + b n sin 6n(p), but the terms of Eq. (S27) are not capable of producing 
so rapidly oscillating functions to linear order. The pinning potential of the spin crystal state is 
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therefore produced by next order terms (sixth order in spin-orbit coupling), for example by 

AF^AsE^KI 2 - (S28) 

q 

We have analyzed the pinning of the spin crystal in the presence of AF^. We obtain for A 3 > 
that the A-crystal gets oriented in such a way that at least two of the six Bragg spots point along 
a (110) direction if the field direction allows for it (i.e., if the field is orthogonal to one of the 
(110) directions). This is compatible with the observations shown in Fig. 2 of the main text. 
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A 



B 




Figure SI: (A) Neutron scattering set-up used in all previous studies to explore the precise 
magnetic structure as function of magnetic field. Note that B is perpendicular to the incident 
neutron beam. (B) Neutron scattering set-up used in our study; the applied magnetic field B 
was parallel to the incident neutron beam. 
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Figure S2: Sketch of the deflection a of the orientation of the skyrmion lattice with respect to 
the sample orientation and magnetic field as observed in sample 1 . This deviation is in perfect 
agreement with the effects of demagnetizing fields. 
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rocking angle (°) 

Figure S3: Typical variation of the scattering intensity as recorded in rocking scans of sample 
2. Intensity represents integrated counts over the individual spots, were labels refer to Fig. 2 in 
the main text. Lines serve as to guide the eye. (A) Intensity in the zero- field-cooled state at 
T = 16 K and B = 0. A spontaneous difference of intensity between the domain populations 
(1,3) and (2,4) is observed. The half- width of the rocking angle implies a magnetic mosaicity 
r] m ps 3.5° consistent with previous work (S32, S3). (B) Intensity in the A-phase on a linear 
scale for spots (5), (6), (7) and (8). The half-width of the rocking angle implies a magnetic 
mosaicity r] m ps 1.75° when taking into account the Lorentz factor of 15% for the location of 
the spots, corresponding to a coherence length £ ps 5500 A ± 10% characteristic of long range 
order. 
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AG 




B I 

Figure S4: Difference of mean-field free energy between various magnetic states and the conical 
state as a function of magnetic field. The energy difference is plotted in units of 7t 2 /4 in the 
limit of small t < 0. 
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UcB, : 



Figure S5: Variation of magnetization amplitudes as a function of field, in the limit of small 
t < 0. The gray region shows the range of |M(r) | inside the unit cell of the A-crystal, plotted in 
units of (— 1/2) 1 / 2 . The two crossing lines correspond to the amplitudes at two specific points in 
the crystal. For the decreasing line, this is the center of an anti-skyrmion (where M points in the 
opposite field direction) and for the increasing line the mid-point between three anti-skyrmions 
(where M points along the field). The dashed line shows the amplitude of the moment in the 
conical phase, which is a constant for all fields. Note that the A-crystal is most favorable, when 
the variation of amplitudes is smallest rather than largest. 



26 




Figure S6: The three contributions to AG = AG mf + AGf« + AGj° § t , plotted as a function 
of 1/A for t = —4 and various fields. The energies per unit volume have been rescaled by 7 or 
A as indicated in the figures. Numerical calculations have been performed up to A = 5. The 
data can be readily extrapolated to the physical cutoff A « 40 (in units of Q). Similar results 
are obtained for smaller values of \t\. For much larger \t\, one would have to increase the cutoff 
to obtain convergence, which makes the numerical computations more demanding. 
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